  subroutine mopy2mo_rhf(icore, itmp, nmo, l, irrepl, isymmo, & 
                   isympopa, isympopb, isymvrta, isymvrtb, nomax, nvmax)   
  use mod_orbit
  use mod_symm
  use mod_iop
  implicit none 
  integer, intent(in) :: nmo, l, irrepl, nvmax, nomax
  integer, intent(in) :: isympopa(nomax, 8), isympopb(nomax, 8)
  integer, intent(in) :: isymvrta(nvmax, 8), isymvrtb(nvmax, 8), isymmo(nmo, 8)
  real*8, intent(in) :: icore(nmo, nmo, nmo)
  real*8, intent(out) :: itmp(*)

  integer i, j, k, a, b, c, d, numka, ioffka 
  integer idx, idxi, idxj, idxk,idxold, idxa, idxb, idxc
  integer noa2, nob2, noab, nva2, nvb2, nvab
  integer nvoa, nvob, nvaob, nvboa, isymabc
  integer nvvabsym(8), nvoabsym(8)
  integer nvobasym(8), nooabsym(8)
  integer nvoaasym(8), nvobbsym(8)
  integer nvv1(8), nvv2(8)
  integer nfirst, nread, nsize, isymoffk(8)
  integer irrep, isym, irrepi, irrepj, irrepij, irrepk
  integer irrepab, irrepa, irrepb, irrepbi, irrepc, irrepaj
  integer irrepci, irrepka, irrepd
  integer ioffij, ioffk, ioffi, ioffab, ioffa, ioffbi
  integer ioffj, ioffc, ioffd, ioffab1, ioffab2
  integer numi, numj, numk, numij, numab, numa, numb
  integer numbi, numc, numab1, numab2, ioffc1, ioffc2
  integer,external :: irpdso, idsymsoc, isymoffso
  real*8, external :: nrm2so, dnrm2
   isym = dirprd(irrepl, 1) 
!  write(6,*)'nirrep', nirrep,l, irrepl
  do irrep = 1, nirrep
     nvvabsym(irrep) = irpdso(irrep, vrta, vrtb, 1)
     nvoabsym(irrep) = irpdso(irrep, vrta, popb, 1)
     nvobasym(irrep) = irpdso(irrep, vrtb, popa, 1)
     nvoaasym(irrep) = irpdso(irrep, vrta, popa, 1)
     nvobbsym(irrep) = irpdso(irrep, vrtb, popb, 1)
     nooabsym(irrep) = irpdso(irrep, popa, popb, 1)
     nvv1(irrep) = irpdso(irrep, vrta, vrta, 0)
     nvv2(irrep) = irpdso(irrep, vrta, vrta, 2)
  enddo
  !ABAB
    !IJKL_ABAB
    idx = 0
    if(l<=popb(irrepl)) then 
!      write(6,*)'l', l, irrepl, isymmo(l, irrepl)
      do irrepk  = 1, nirrep 
         irrepij = dirprd(irrepk, isym) 
         numk = popa(irrepk)
         ioffk = isymoffso(irrepk, isym, nooabsym, popa, 1)
         numij = irpdso(irrepij, popa, popb, 1)
         do k = 1, numk 
         do irrepj = 1, nirrep
            irrepi = dirprd(irrepj, irrepij) 
            numj = popb(irrepj) 
            numi = popa(irrepi)
            ioffij = isymoffso(irrepj, irrepij, popa, popb, 1)
            do j = 1, numj 
            do i = 1, numi
               idxi = isympopa(i, irrepi) 
               idxj = isympopb(j, irrepj) 
               idxk = isympopa(k, irrepk) 
               !<Ij||Kl>=(IK|jk)
               idxold = idx 
               idx = (ioffk-1)+ (k-1)*numij + ioffij-1 + (j-1)*numi + i    
               itmp(idx) = icore(idxi,idxk,idxj) 
!              if(idx-idxold/=1) then 
!                 write(6,'(A, 6I6)')'error in idx of ijkl', ioffk, isymmo(l, irrepl), idx, idxold
!              endif
            enddo
            enddo
            enddo
         enddo
      enddo
      nsize = idsymsoc(isym, popa, popb, 1, popa, zorb, 1)
      nsize = irpdso(isym, nooabsym, popa, 1) 
!      write(6,'(A, I5, I5, I5, F18.8)')'norm of ijkl',  isymmo(l, irrepl), isym, irrepl, nrm2so(nsize, itmp, 1) 
      call putpqrso(itmp, l, irrepl, 1,  'IJKLAB__')
    endif
!    write(6,*)'finished IJKLAB' 
    !ABIJ_ABAB
    if(l<=popb(irrepl)) then
      j = l 
      irrepj = irrepl 
      do irrepi = 1, nirrep
         irrepab = dirprd(irrepi, isym) 
         numi = popa(irrepi) 
         numab = irpdso(irrepab, vrta, vrtb, 1)
         ioffi = isymoffso(irrepi, isym, nvvabsym, popa, 1) 
         do irrepb = 1, nirrep
            irrepa = dirprd(irrepb, irrepab)
            numb = vrtb(irrepb) 
            numa = vrta(irrepa)
            ioffab = isymoffso(irrepb, irrepab, vrta, vrtb, 1)
            do i = 1, numi
            do b = 1, numb 
            do a = 1, numa
               idxi = isympopa(i, irrepi) 
               idxa = isymvrta(a, irrepa) 
               idxb = isymvrtb(b, irrepb) 
               idx = (ioffi-1) + (i-1)*numab + ioffab-1 + (b-1)*numa + a
               itmp(idx) = icore(idxa, idxi, idxb)
            enddo
            enddo
            enddo
          enddo
       enddo
!       write(6,*)'l,irrepl', l, irrepl
       call putpqrso(itmp, l, irrepl, 1,  'ABIJAB__')
    endif
!    write(6,*)'finished ABIJAB' 
!   IABJ_AAAA = BIAJ
! <IA||BJ>=(IB|AJ)-(AB|IJ)
    if(l<=popb(irrepl)) then 
       j = l; irrepj = irrepl 
      do irrepa = 1, nirrep
         irrepbi = dirprd(irrepa, isym) 
         numa = vrtb(irrepa) 
         numbi = irpdso(irrepbi, vrta, popa, 1) 
         ioffa = isymoffso(irrepa, isym, nvoaasym, vrta, 1)
         do irrepi = 1, nirrep 
            irrepb = dirprd(irrepi, irrepbi) 
            numi = popa(irrepi) 
            numb = vrta(irrepb) 
            ioffbi = isymoffso(irrepi, irrepbi, vrta, popa, 1)
            do a = 1, numa
            do i = 1, numi
            do b = 1, numb
               idxa = isymvrta(a, irrepa)
               idxi = isympopa(i, irrepi)
               idxb = isymvrta(b, irrepb)
               idx = (ioffa-1) + (a-1)*numbi + ioffbi - 1 +(i-1)*numb+b
               itmp(idx) = icore(idxi, idxb, idxa) - icore(idxa, idxb, idxi)
            enddo
            enddo
            enddo
          enddo
        enddo
       call putpqrso(itmp, l, irrepl, 1,  'BIAJAAAA')
      endif
!    write(6,*)'finished BIAJ1'

!   !IaBj_ABAB = BIaj
    if(l<=popb(irrepl)) then 
       j = l; irrepj = irrepl 
      do irrepa = 1, nirrep
         irrepbi = dirprd(irrepa, isym) 
         numa = vrtb(irrepa) 
         numbi = irpdso(irrepbi, vrta, popa, 1) 
         ioffa = isymoffso(irrepa, isym, nvoaasym, vrtb, 1)
         do irrepi = 1, nirrep 
            irrepb = dirprd(irrepi, irrepbi) 
            numi = popa(irrepi) 
            numb = vrta(irrepb) 
            ioffbi = isymoffso(irrepi, irrepbi, vrta, popa, 1) 
            do a = 1, numa 
            do i = 1, numi 
            do b = 1, numb 
               idxa = isymvrtb(a, irrepa) 
               idxi = isympopa(i, irrepi) 
               idxb = isymvrta(b, irrepb) 
               idx = (ioffa-1) + (a-1)*numbi + ioffbi - 1 +(i-1)*numb+b
               itmp(idx) = icore(idxi, idxb, idxa) 
            enddo 
            enddo
            enddo
          enddo
        enddo
       call putpqrso(itmp, l, irrepl, 1,  'BIAJAABB')
      endif
!    write(6,*)'finished BIAJ2'

!   !iABj_BAAB = BiAj
    if(l<=popb(irrepl)) then
       j = l; irrepj = irrepl
      do irrepa = 1, nirrep
         irrepbi = dirprd(irrepa, isym)
         numa = vrta(irrepa)
         numbi = irpdso(irrepbi, vrta, popb, 1)
         ioffa = isymoffso(irrepa, isym, nvoaasym, vrta, 1)
         do irrepi = 1, nirrep
            irrepb = dirprd(irrepi, irrepbi)
            numi = popb(irrepi)
            numb = vrta(irrepb)
            ioffbi = isymoffso(irrepi, irrepbi, vrta, popb, 1)
            do a = 1, numa
            do i = 1, numi
            do b = 1, numb
               idxa = isymvrta(a, irrepa)
               idxi = isympopb(i, irrepi)
               idxb = isymvrta(b, irrepb)
               idx = (ioffa-1) + (a-1)*numbi + ioffbi - 1 +(i-1)*numb+b
               itmp(idx) = -icore(idxa, idxb, idxi)
            enddo
            enddo
            enddo
          enddo
        enddo
       call putpqrso(itmp, l, irrepl, 1,  'BIAJABAB')
      endif

!    write(6,*)'finished BIAJ3'
!   !KAIJ_ABAB
    if(l<=popb(irrepl)) then
      j = l
      irrepj = irrepl
      do irrepi = 1, nirrep
         irrepka = dirprd(irrepi, isym)
         numi = popa(irrepi)
         numka = irpdso(irrepka, popa, vrtb, 1)
         ioffi = isymoffso(irrepi, isym, nvoabsym, popa, 1)
         do irrepa = 1, nirrep
            irrepk = dirprd(irrepa, irrepka)
            numa = vrtb(irrepa)
            numk = popa(irrepk)
            ioffka = isymoffso(irrepa, irrepka, popa, vrtb, 1)
            do i = 1, numi
            do a = 1, numa
            do k = 1, numk
               idxi = isympopa(i, irrepi)
               idxa = isymvrta(a, irrepa)
               idxk = isympopb(k, irrepk)
               idx = (ioffi-1) + (i-1)*numka + ioffka-1 + (a-1)*numk + k
               itmp(idx) = icore(idxk, idxi, idxa)
            enddo
            enddo
            enddo
          enddo
       enddo
!       write(6,*)'l,irrepl', l, irrepl
        call putpqrso(itmp, l, irrepl, 1,  'KAIJAB__')
    endif

    if(eri_level==0.or.eri_level==1) then 
!   !ABCI_ABAB
      if(l<=popb(irrepl)) then 
      do irrepc = 1, nirrep
         irrepab = dirprd(irrepc, isym) 
         numc = vrta(irrepc) 
         numab = irpdso(irrepab, vrta, vrtb, 1) 
         ioffc = isymoffso(irrepc, isym, nvvabsym, vrta, 1) 
         do irrepb = 1, nirrep
            irrepa = dirprd(irrepb, irrepab) 
            ioffab = isymoffso(irrepb, irrepab, vrta, vrtb, 1) 
            numa = vrta(irrepa)
            numb = vrtb(irrepb) 
            do c = 1, numc
            do b = 1, numb
            do a = 1, numa   
               idxc = isymvrta(c, irrepc) 
               idxb = isymvrtb(b, irrepb) 
               idxa = isymvrta(a, irrepa) 
               idx = ioffc - 1 + (c-1)*numab + ioffab - 1 + (b-1)*numa + a 
               itmp(idx) = icore(idxa, idxc, idxb) 
             enddo
             enddo
             enddo
           enddo
        enddo
        call putpqrso(itmp, l, irrepl, 1,  'ABCIAB__')
       endif 
     endif
!    write(6,*)'finished ABCI'

!   !IJKA_ABAB
!   if(l>popb(irrepl)) then
!     a = l - popb(irrepl) 
!     irrepa = irrepl
!     do irrepk  = 1, nirrep
!        irrepij = dirprd(irrepk, isym)
!        numk = popa(irrepk)
!        ioffk = isymoffso(irrepk, isym, nooabsym, popa, 1)
!        numij = irpdso(irrepij, popa, popb, 1)
!        do k = 1, numk
!        do irrepj = 1, nirrep
!           irrepi = dirprd(irrepj, irrepij)
!           numj = popb(irrepj)
!           numi = popa(irrepi)
!           ioffij = isymoffso(irrepj, irrepij, popa, popb, 1)
!           do j = 1, numj
!           do i = 1, numi
!              idxi = isympopa(i, irrepi)
!              idxj = isympopb(j, irrepj)
!              idxk = isympopa(k, irrepk)
!              !<Ij||Ka>=(IK|ja)
!              idxold = idx
!              idx = (ioffk-1)+ (k-1)*numij + ioffij-1 + (j-1)*numi + i
!              itmp(idx) = icore(idxi,idxk,idxj)
!           enddo
!           enddo
!           enddo
!        enddo
!     enddo
!      write(6,*)'a, irrepa', a, irrepa
!     call putpqrso(itmp, a, irrepa, 1,  'IJKAAB__')
!   endif
!    write(6,*)'finished IJKA'


!   !ABCD1 and ABCD2
!   !<AB||CD> = (AC|BD)-(BC|AD)
!   ABCD2 = (AC|BD)+(BC|AD)
    if(eri_cfour==1) then
!   !<Ab||Cd> = (AC|bd)
      if(l>popb(irrepl)) then
      d = l - popb(irrepl)
      irrepd = irrepl
      do irrepc = 1, nirrep
         irrepab = dirprd(irrepc, isym)
         numc = vrta(irrepc)
         numab = irpdso(irrepab, vrta, vrtb, 1)
         ioffc = isymoffso(irrepc, isym, nvvabsym, vrta, 1)
         do irrepb = 1, nirrep
            irrepa = dirprd(irrepb, irrepab)
            ioffab = isymoffso(irrepb, irrepab, vrta, vrtb, 1)
            numa = vrta(irrepa)
            numb = vrtb(irrepb)
            do c = 1, numc
            do b = 1, numb
            do a = 1, numa
               idxc = isymvrta(c, irrepc)
               idxb = isymvrtb(b, irrepb)
               idxa = isymvrta(a, irrepa)
               idx = ioffc - 1 + (c-1)*numab + ioffab - 1 + (b-1)*numa + a
               itmp(idx) = icore(idxa, idxc, idxb)
             enddo
             enddo
             enddo
           enddo
        enddo
        call putpqrso(itmp, d, irrepd, 1,  'ABCDAB__')
       endif
    endif
 
    if(eri_level==0) then
       if(l>popb(irrepl)) then
       d = l - popb(irrepl)
       irrepd = irrepl
       do irrepc = 1, nirrep
          irrepab = dirprd(irrepc, isym)
          numc = vrta(irrepc)
          numab = irpdso(irrepab, vrta, vrta, 0)
          do irrepb = 1, nirrep
             irrepa = dirprd(irrepb, irrepab)
             ioffc = isymoffso(irrepc, isym, nvv1, vrta, 1)
             ioffab = isymoffso(irrepb, irrepab, vrta, vrta, 0)
             numa = vrta(irrepa)
             numb = vrtb(irrepb)
             if(irrepa>irrepb) cycle 
             if(irrepa==irrepb) then 
               do c = 1, numc
               do b = 2, numb
               do a = 1, b-1
                  idxc = isymvrta(c, irrepc)
                  idxb = isymvrtb(b, irrepb)
                  idxa = isymvrta(a, irrepa)
                  idx = ioffc - 1 + (c-1)*numab + ioffab - 1 + (b-1)*(b-2)/2 + a
                  itmp(idx) = icore(idxa, idxc, idxb) - icore(idxb, idxc, idxa)
                enddo
                enddo
                enddo
             elseif(irrepa<irrepb) then 
               do c = 1, numc
               do b = 1, numb
               do a = 1, numa 
                  idxc = isymvrta(c, irrepc)
                  idxb = isymvrtb(b, irrepb)
                  idxa = isymvrta(a, irrepa)
                  idx = ioffc - 1 + (c-1)*numab + ioffab - 1 + (b-1)*numa + a
                  itmp(idx) = icore(idxa, idxc, idxb) - icore(idxb, idxc, idxa)
                enddo
                enddo
                enddo
              endif
            enddo
         enddo
         call putpqrso(itmp, d, irrepd, 1,  'ABCDAA__')
        endif
!!!!    write(6,*)'finished abcdaa'
       if(l>popb(irrepl)) then
       d = l - popb(irrepl)
       irrepd = irrepl
       do irrepc = 1, nirrep
          irrepab = dirprd(irrepc, isym)
          numc = vrta(irrepc)
          numab = irpdso(irrepab, vrta, vrta, 2)
          do irrepb = 1, nirrep
             irrepa = dirprd(irrepb, irrepab)
             ioffc = isymoffso(irrepc, isym, nvv2, vrta, 1)
             ioffab = isymoffso(irrepb, irrepab, vrta, vrta, 2)
             numa = vrta(irrepa)
             numb = vrtb(irrepb)
             if(irrepa>irrepb) cycle
             if(irrepa==irrepb) then
               do c = 1, numc
               do b = 1, numb
               do a = 1, b
                  idxc = isymvrta(c, irrepc)
                  idxb = isymvrtb(b, irrepb)
                  idxa = isymvrta(a, irrepa)
                  idx = ioffc - 1 + (c-1)*numab + ioffab - 1 + (b-1)*b/2 + a
                  itmp(idx) = icore(idxa, idxc, idxb) + icore(idxb, idxc, idxa)
                enddo
                enddo
                enddo
             elseif(irrepa<irrepb) then
               do c = 1, numc
               do b = 1, numb
               do a = 1, numa
                  idxc = isymvrta(c, irrepc)
                  idxb = isymvrtb(b, irrepb)
                  idxa = isymvrta(a, irrepa)
                  idx = ioffc - 1 + (c-1)*numab + ioffab - 1 + (b-1)*numa + a
                  itmp(idx) = icore(idxa, idxc, idxb) + icore(idxb, idxc, idxa)
                enddo
                enddo
                enddo
              endif
            enddo
         enddo
         do irrepc = 1, nirrep 
            irrepab = dirprd(irrepc, isym) 
            ioffc = isymoffso(irrepc, isym, nvv2, vrta, 1) 
            numc = vrta(irrepc) 
            numab = irpdso(irrepab, vrta, vrta, 2) 
            if(irrepd==irrepc) then 
              do c = 1, numc 
                 if(c==d) then 
                   ioffab=ioffc + (c-1)*numab 
                   call scalso(numab, 0.5d0, itmp(ioffab), 1)
                 endif
              enddo
            endif
          enddo
         call putpqrso(itmp, d, irrepd, 1,  'ABCDAB2_')
        endif
   endif

  return
  end
    
  
